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We treat three-dimensional bosonic clusters with up to TV — 40 atoms, interacting additively 
through two-body Van der Waals potentials, in the near-threshold regime. Our study includes 
super-borromean systems with TV atoms for which all subsystems are unbound. We determine the 
energetics and structural properties such as the expectation value of the interparticle distance as 
a function of the coupling strength. It has been shown that the coupling strength gi , for which 
the TV-body system becomes unbound, is bounded by the coupling constant gi N , for which the 



> (TV - l)gi N ~ 1] /N. By fitting 



next smaller system with TV — 1 atoms becomes unbound, i.e., g 
our numerically determined ground state energies to a simple functional form with three fitting 
parameters, we determine the relationship between g[ N ^ and g* _1 \ Our trimer and tetramer 
energies fall on the so-called Tjon line, which has been studied in nuclear physics. We confirm 
the existence of generalized Tjon lines for larger clusters. Signatures of the universal behavior of 
weakly-bound three-dimensional clusters can possibly be observed in ultracold Bose gases. 
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I. INTRODUCTION 

Weakly-bound few-body systems have been studied ex- 
tensively by the atomic, nuclear and condensed matter 
physics community since the early days of quantum me- 
chanics. Within the framework of non-relativistic quan- 
tum mechanics, the properties of a many-body system are 
determined by the mass and statistics of the constituents 
and by the potential energy surface. In many cases, the 
many-body potential energy surface can be approximated 
quite accurately by a sum of two-body potentials. The 
interaction strength g of the TV-body system is then de- 
termined by the underlying two-body potential and, as- 
suming TV identical mass m particles, the mass m. The 
critical coupling constant g* N \ for which the TV-body 
cluster becomes unbound, defines the threshold. This 
paper investigates the near-threshold regime of bosonic 
three-dimensional TV-body clusters, i.e., the regime where 
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> 



g{ N \ This near-threshold regime is particularly 
interesting since some properties of the bosonic many- 
body system become independent of the details of the 
underlying potential energy surface, i.e., some properties 
of weakly-bound clusters consisting of TV bosons become 
universal as g -> gi N) 0, S S S S E II II @ 

Our study includes the characterization of "super- 
borromean" TV-body clusters 01 ■ Borromean trimers, 
which consist of three bosons for which each dimeric 
subsystem is unbound, have been studied in detail in 
the literature 0, 0, Super-borromean clusters, 

which consist of TV bosons for which all subsystems with 
TV — n, where n = 1, • • • , N — 2, are unbound, in con- 
trast, have not been studied in much detail. To char- 
acterize these delicate systems, we perform precise dif- 
fusion Monte Carlo (DMC) calculations for clusters in- 
teracting additively through realistic shape-dependent 
two-body Van der Waals potentials. We determine the 



critical coupling strengths gi N ^ for atomic clusters with 
up to TV = 40 bosons and compare with variational 
bounds. The near-threshold behavior of weakly-bound 
three-dimensional bosonic clusters has been investigated 
in a series of papers before We believe, how- 

ever, that advances in the theoretical understanding, in- 
cluding predictions derived using effective theories and 
zero-range models H, and in the numerical treat- 
ment make it worthwhile to revisit the characterization 
of weakly-bound three-dimensional clusters. In particu- 
lar, we present more accurate energies for a larger range 
of coupling strengths and for a larger range of cluster 
sizes than previous studies. 

The present study is additionally motivated by recent 
experiments on extremely weakly-bound molecules cre- 
ated from ultracold Bose and Fermi gases. Utilizing 
Feshbach resonances the effective interaction strength 
between two atoms at ultracold temperatures can be 
changed essentially at will through application of an ex- 
ternal magnetic field [lol ITi| . The existence of this exter- 
nal knob has led to the observation of extremely weakly- 
bound diatomic molecules in highly-excited vibrational 
states [13, HM llil an d provided evidence for the forma- 
tion of Efimov trimer states [53, in an ultracold en- 
vironment. Furthermore, recent experiments on cold Cs 
atoms evidence the creation of larger weakly-bound clus- 
ters |22j : these experiments point towards Feshbach en- 
gineering of weakly-bound clusters. Feshbach resonances 
arise from the coupling of two Born-Oppenheimer poten- 
tial curves through a hyperfine Hamiltonian and require, 
in general, a multi-channel description. In the case of 
a broad resonance, however, the change of the effective 
scattering length can be described within a single channel 
model 23] . Using a single channel approximation, this 
paper describes weakly-bound three-dimensional bosonic 
clusters with varying atom-atom scattering lengths with 
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up to N = 40 atoms. 

Section fll Al dcscribes the many-body Hamiltonian and 
the characteristics of the underlying two-body potential. 
Section III Bl is devoted to a discussion of our numerical 
approach to solving the many-body Schrodinger equa- 
tion. Our results for the energetics and structural prop- 
erties are presented in Sees. Hill and IIVI respectively, and 
our conclusions in Sec. 



II. SYSTEM AND NUMERICAL APPROACH 



A. Many-body Hamiltonian 



Consider the Hamiltonian H for N bosons with mass 
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V(r jk ) 



(1) 



where V| and 



rjk denote respectively the 3D Laplace 
operator of the jth boson and the internuclear dis- 
tance between particles j and k. This Hamiltonian 
assumes a many-body potential energy surface written 
as a sum of atom-atom potentials V(r). Our calcula- 
tions are performed for a realistic Van der Waals triplet 
tritium-tritium potential [H El HI IHIH HJ , which 
is repulsive at short interparticle distances r and falls 



off at large r as Xm= 



6,8,— 



-C n r" n . Figure ^ shows 
the tritium-tritium potential as a function of the inter- 
particle distance r. The potential has a minimum of 
depth — 4.6cm -1 at r w 7.8ao, where ao denotes the 
Bohr radius. Solving the one-dimensional scaled radial 
Schrodinger equation shows that the tritium dimer has 
no bound state 
length a, 



Internuclear distance r (in units a ) 



FIG. 1: Triplet tritium-tritium interaction potential as a func- 
tion of the internuclear distance r. 



As alluded to in the introduction, the coupling strength 
can be varied experimentally via a Feshbach reso- 
nance [30l l3l| . Although a full description of a Fesh- 
bach resonance requires the coupling between at least 
two channels — in tritium, e.g., of the singlet and triplet 
potential curves (cou pled through a long-range hyper- 
fine Hamiltonian) |lOl |32| — , some properties can be de- 
scribed within a single channel model. We thus envi- 
sion that changing the atom mass in our single channel 
treatment can be mapped to changing the strength of 
an external magnetic field, and hence of the atom-atom 
scattering length, in the vicinity of a two-body Feshbach 
resonance. We expect that our calculations uncover the 
"generic" behaviors of three-dimensional bosonic Van der 
Waals clusters, which are interacting additively through 
two-body potentials with repulsive short-range core and 



(see also Sec. |HTAj|. The scattering long-range tail with leading —Ce/r 6 term. In particular, 



a = lim — 



tan((5(fc)) 



(2) 



of the tritium dimer is negative, i.e., a ~ -82. la [13, 
which indicates that the dimer is only slightly short of 
binding. In Eq. (J2J), S(k) denotes the energy-dependent 
s-wave phase shift and k the relative wave vector of the 
equivalent one-body problem with reduced mass m/2. 

Our interest in this paper is in a detailed description 
of weakly-bound clusters with varying coupling constant 
near threshold. To change the coupling strength g of the 
cluster, we vary the atom mass to, i.e., we consider "ar- 
tificial" clusters with atom masses that are heavier and 
lighter than the tritium mass. By rewriting the many- 
body Schrodinger equation in scaled units, it can be read- 
ily seen that changing the atom mass changes the cou- 
pling strength. For example, for systems interacting ad- 
ditively through Lenard-Jones potentials with well depth 
e and length scale a the coupling constant g is directly 
proportional to the atom mass, g = Ameer 2 /h 2 . 



we believe that usage of a different two-body potential 
in Eq. JQl will result in the same qualitative but possibly 
different quantitative behaviors of weakly-bound bosonic 
clusters. 

Pluses in Fig. [21 show the atom- atom scattering length 
a for the tritium-tritium potential as a function of the 
atom mass to. The scattering length a diverges at 
to w 5933.4(2)m e , where m e denotes the electron mass 
and the value in round brackets the uncertainty of m 
arising from the numerical determination of the scatter- 
ing length a. Since this is the mass at which the dimer 
becomes unbound, we refer to this mass as the critical 
mass of the dimer. A vertical solid line in Fig. [21 
marks the value of to; . We find that the scattering 
length a vanishes for m ps 2311.0(2)TO e . This mass value 
is indicated by a vertical solid line in Fig. [21 and puts 
an upper bound on the critical mass for the bulk system 
(N -> 00), i.e., to! oo) < 2311. 0to £ 1]. This bound is 
obtained variationally by expanding the energy in terms 
of the density. For negative a, the leading order in the 
expansion becomes negative and the bulk system is nec- 
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FIG. 2: Pluses and diamonds show the atom-atom scat- 
tering length a and the effective range r e ff, respectively, 
as a function of the atom mass m for the triplet tritium- 
tritium potential (see text). To guide the eye dotted lines 
connect the symbols. Vertical solid lines indicate the mass 
m ~ 2311. 0m e , at which a goes through zero, and the criti- 

(2) 

cal mass ml ~ 5933.4m,,, at which the scattering length a 
diverges. 



essarily bound |l| . Section IIII CI compares our critical 
masses mi N ^ calculated for up to N — 40 atoms with the 
variational upper bound for ml 00 -* . 

Diamonds in Fig.[21show the effective range J~ e //, which 
we calculate through the relationship 



a(k) 



(3) 



as a function of the atom mass to. In Eq. a(k) 
denotes the energy-dependent scattering length, a(k) = 
— tan(<5(fc))/fc. The effective range is largest in the re- 
gion where the two-body scattering length vanishes. In 
the region where a diverges, r e ff takes on values of the 
order of lOao- The scattering length a, the effective 
range r e // and the Van der Waals length rydw, where 
rvdw = (mCe/H 2 ) 1 / 4 , are the relevant length scales of 
the two-body problem near threshold. For a zero-range 
model with a single parameter, namely the scattering 
length a, to be applicable for N = 2, a needs to be the 
largest length scale in the problem. This condition can 
be expressed as 



\E 2 \ <C min 



(4) 



where E 2 denotes the ground state energy of the dimer. 
For zero-range models to be applicable to clusters with 
N > 2, a condition similar to Eq. Q, possibly with an 
additional scaling factor N or N(N — l)/2, needs to be 
fulfilled. 



B. Numerical treatment of iV-body clusters 

To determine the ground state energy and wave func- 
tion of the two-body system [Eq. with N = 2], we 
separate off the center of mass motion and scale the 
wave function for the interparticle distance to remove 
first derivative terms in the kinetic energy operator. The 
scaled one-dimensional radial Schrodinger equation can 
then be solved by diagonalizing the Hamiltonian using 
B-splines. To treat very weakly-bound dimers with vary- 
ing mass to, we optimize the adaptive grid (i.e., the grid 
spacings, the number of grid points and the integration 
interval) for each mass. The upper integration limit is 
determined by the size of the bound state; integrating 
the Schrodinger equation out to roughly 100a leads to 
converged results for all two-body systems considered in 

sec, min 

The calculations of the trimer energies are, due to the 
larger number of degrees of freedom, more involved than 
those of the dimer energies. Separating off the center of 
mass motion reduces the nine-dimensional problem to a 
six-dimensional problem. Since we are in this paper pri- 
marily interested in ground state properties, we restrict 
ourselves to states with vanishing total angular momen- 
tum, i.e., J = 0. The resulting three-dimensional Hamil- 
tonian can be written in terms of Whitten and Smith' 
hyperspherical coordinates [3j|, which allow the Bosc 
symmetry to be accounted for readily. To solve the cor- 
responding scaled Schrodinger equation, we expand the 
wave function in angle-dependent channel functions $, 
which depend parametrically on the hyperradius R, and 
a set of weight functions F(R) . Our numerical implcmcn- 
tion is described in Ref. j34j. We check the convergence 
by changing the hyperangular grid, the hyperradial grid, 
the number of channel functions included in the expan- 
sion, and the step size used in the numerical determi- 
nation of the derivatives of the channel functions. The 
trimer ground state energies presented in Sec. IIII 51 have 
an accuracy of a few percent. For selected trimers, we 
also report the first excited state energy with J = 0. 

For cluster systems with more than a few atoms, 
basis set expansion-type techniques become computa- 
tionally unfeasible. Consequently, we solve the many- 
body Schrodinger equation for TV > 4 using alternative 
techniques, i.e., the variational quantum Monte Carlo 
(VMC) method and the DMC method with importance 
samp ling |35j ]. Our numerical implementations follow 
Ref. [3g. The VMC method minimizes the energy of 
the cluster system by optimizing the many-body wave 
function, which is written in terms of a set of parame- 
ters p. The optimized variational wave function ipx then 
enters our DMC calculations, which result in essentially 
exact ground state energies, as a guiding function. Due 
to the stochastic nature of the MC algorithms, the DMC 
energies reported in Sec. IIII Ol have statistical uncertain- 
ties. 

We use two different functional forms for the varia- 
tional wave function. For small clusters with about up 
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to N — 10 atoms, each atom has roughly the same av- 
erage distance to all other atoms in the cluster. In this 
case, our variational wave function ipT is written as a 
product of pair wave functions cj> |37| , 



,r N ) 
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where 



4>(r) = exp 



P-5 P2 , , , 

— -Polog(r) 



Pir 



(5) 



(6) 



For larger clusters, the variational wave function given in 
Eqs. © and © does not give a good variational energy 
and we additionally include a variational Fermi function 
which depends on the distance Rj of the jth atom to the 
center of mass of the cluster |3g , 



N 

n< 



(rjk) 



N 



(7) 



Here, (f> is given by Eq. JBJ with pg = p\ = and 

-l 



4>{R) 



1 + exp 



R-Pe 
Pa 



(8) 



The variational parameters p e and p a determine the size 
of the cluster and the sharpness of the cluster's surface 
region, respectively. For each cluster system considered, 
we optimize the variational parameters by minimizing the 
energy expectation value {ipr\H \4>t) / ('4>t\'4't} ■ The 3N- 
dimensional integrals are evaluated using the Metropolis 
algorithm. Our VMC energies, except for those for the 
systems closest to threshold (see Sec. lIIlCj) . recover more 
than about 75-80% of the essentially exact DMC ground 
state energies. 

The DMC calculations become computationally more 
demanding as we approach the threshold, since the ki- 
netic and potential energy nearly cancel. In the near- 
threshold regime great care has to be taken to avoid any 
guiding function bias and to ensure convergence of the 
DMC calculations. To check that our DMC code de- 
scribes extremely weakly-bound clusters accurately, we 
compare the DMC energies for the trimer with those cal- 
culated by the hyperspherical B-splinc treatment. We 
find agreement to within the statistical uncertainty for 
m > 6000m e but do not obtain reliable DMC energies 
for significantly smaller masses. 

Unlike the DMC energy expectation value, which is 
essentially exact (except for statistical uncertainties and 
possible time step errors), the expectation value of any 
structural quantity B is in the "standard" DMC algo- 
rithm calculated with respect to the mixed density, 



{B) DMC = WSIV'T)/ Wr). 



(9) 



Here, denotes the exact stationary ground state wave 
function |35|. To improve upon this mixed estimator, 



we calculate the so-called extrapolated expectation value 

{B) ex m, 

(B) ex = 2{B) DMC - (B)vmc, (10) 
where {B)vmc denotes the VMC expectation value, 

(B)vmc = (lMB|lM/GMVr>- ( n ) 

For the systems studied in this paper, the extrapolated 
expectation values {B) ex are expected to be fairly close 
to the exact expectation values. Section ITVl reports ex- 
pectation values for the pair distribution function P(r) 
and the interparticle distance r. 

III. ENERGETICS 

This section presents our numerically determined en- 
ergies for clusters with up to 40 atoms and their inter- 
pretation. 

A. N = 2 

Pluses in Fig. |3{a) show the absolute value of the s- 
wave ground state energies Ei for two particles inter- 
acting through the triplet tritium-tritium potential for 
a number of different m, i.e., m G [5933. 4m e , IOOOOtoJ. 
The tritium dimer itself is, as mentioned in Sec. Ill Al 
unbound. The ground state energies shown in Fig. 
extend over nearly ten orders of magnitude; Ei for the 
most weakly-bound dimer considered with m — 5933.4 is 
—6 x 10~ 11 cm~ 11 , and that for the most strongly-bound 
dimer considered with m = 10000m e is — 0.19cm -1 . 

We can compare the numerically determined ground 
state energies E 2 with the energies E^ predicted from 
a zero-range model, which supports a bound state for 
positive a, 



E' 



(12) 



A solid line in Fig. Ota) shows the absolute value of the 
zero-range energies E\. In the region where the zero- 
range model provides a good description of the dimer 
energies, the scattering length is the largest length scale 
in the problem, i.e., a r e ff and a ryaw (see Fig- EJ • 
As a becomes comparable to r e f f , the energies E\ deviate 
visibly from the exact energies E 2 and Eq. (|12|l has to be 
modified to account for the dependence of the energy on 
the effective range r e // in addition to a. 

We now describe an analysis that allows an accurate 
determination of the critical mass ml from the two- 
body ground state energies. In principle, this analysis is 
not needed since our scattering length calculations allow 

(2) 

the critical mass ml ' to be determined with high accu- 
racy (see Sec. Ill A)) . The analysis presented in the next 
two paragraphs for N = 2 is meant as a test of principle; 
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FIG. 3: (a) Pluses show the absolute value of the numerically 
determined two-body ground state energies E2 as a function 
of the atom mass m. To guide the eye, a dotted line connects 
the symbols. For comparison, a solid line shows the absolute 
value of the two-body energies E\ , Eq. 1121 , obtained from the 
scattering length a. (b) Pluses show the scaled dimer ground 
state energies y/ m\Ei\ as a function of the atom mass m. 
The scaled energies \J m\Ei\ vary to first order linearly with 
m. A solid line shows our fit to the scaled energies y/m\E2\, 
treating c^ 2 ', and ml 2 ' as fitting parameters (see text), 
(c) Solid, dashed and dotted lines show the lowest two-body 
energies for / = 0, 1 and 2, respectively, as a function of the 
atom mass m. Note that the mass range shown in the lowest 
panel differs from the mass ranges shown in the upper two 
panels. 



an analogous analysis is in Sees. IIIIBl and ITU CI applied 
to larger clusters. For N — 3, accurate calculations of 
the dimer plus atom scattering length can be performed 
but are not pursued here. For N > 3, however, only ap- 
proximate calculations for the N — 1 plus atom scattering 
length have been performed to date [5J, |4fJ . 

Within effective range theory, the two-body ground 
state energies E2 near threshold are determined by 

Taylor-expanding the logarithmic derivative of the bound 

(2 

state wave function about the critical mass ml 



v ^4^| = ^cf ) (m-ml 2) ) i 

i=l 

,(2) 



(13) 



N 


ml [m e \ 


c i [\Zcm- 1 /m e ] 


4 [y/ dor 1 /ml] 


2 


5933.5(1) 


0.011623(2) 


-2.243(5) x 10~ 7 


3 


5352(17) 


0.026(2) 


-2.6(6) x 10~ 6 


4 


4836(08) 


0.033(1) 


-1.9(2) x 10~ 6 


5 


4527(15) 


0.042(2) 


-3.2(6) x 10~ 6 


6 


4311(18) 


0.051(2) 


-4.2(9) x 10~ 6 


7 


4097(06) 


0.056(1) 


-3.7(4) x 10~ 6 


8 


3992(09) 


0.067(2) 


-7.0(9) x 10~ 6 


9 


3867(16) 


0.072(3) 


-6.6(1.7) x 10~ 6 


10 


3789(14) 


0.084(4) 


-1.0(3) x 10~ 5 


20 


3142(26) 


0.090(8) 


-1.4(6) x 10~ 5 


40 


2919(13) 


0.202(7) 


-2.9(6) x 10" 5 



TABLE I: Fitting parameters m* , c[ 
dimensional bosonic clusters with N = 2 — 10, 20 and 40 inter- 
acting additively through a triplet tritium-tritium potential. 
The numbers in brackets indicate the uncertainties of the fit, 
neglecting possible uncertainties of the energies (see text). 



(JV) 

=1 



J*0 



in Fig. EJb) show the scaled energies y/rn\E^, which 
vary roughly linearly with m. Close inspection, how- 
ever, reveals deviations from a linear behavior. This 
indicates that the first term in the expansion given 

by Eq. (|13fl is dominant, but that the second expan- 

(2) 

sion coefficient c 2 contributes non-ncglegibly. To dc- 



,(2) 



termine ml"' , c\" and cv, , we fit our scaled energies 
for m £ [5933. 4m e , lOOOOmJ to the first two terms of 
the expansion given by Eq. l(T3|). The resulting fit with 

to1 2) = 5933.5(l)m e , cf } = 1.1623(2) x 10~ 2 ^cm" 1 Jm e 

and 4 2) = -2.243(5) x 10" 7 Vcm" 1 Jm% (see Table HJ) 
agrees well with the exact energies and is shown by a solid 
line in Fig. OJb). The numbers in round brackets indi- 
cate the uncertainty of the fitting parameters, excluding 
possible numerical inaccuracies of the two-body energies. 
The critical mass extracted by fitting to the dimer bound 

state energies is in excellent agreement with the critical 

(21 

mass ml = 5933.4(2)m e determined from the scattering 
length calculations (see Sec. Ill All . 

The calculations for larger clusters necessarily cover, 
due to numerical difficulties, a smaller range of energies, 
i.e., we are not able to perform bound state calculations 
as close to threshold as for the dimer. Furthermore, our 
cluster energies for N > 3 can only be determined within 
a statistical uncertainty, which adds an additional com- 
plication. If we exclude the two-body energies very close 
to threshold from our fit, i.e., if we perform a fit to the 
scaled energies with m G [6800?7j e , lOOOOmJ (which is 
roughly comparable to the corresponding ranges consid- 
ered for N > 3, see Sees. HTlBl and lilTcl . we find a 



critical mass 



,(2) 



5932.1(8)m e , where the uncertainty 



where the c\ denote expansion coefficients. Pluses 



in brackets reflects, as above, the uncertainty of the fit. 
Since c 2 is negative the critical mass predicted by this 
fit, which excludes the energies closest to threshold, is 
expected to be smaller than the exact threshold value. 
The deviation from the fit that includes the whole mass 
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FIG. 4: Pluses and diamonds show respectively the absolute 
value of the ground state and first excited state energies with 
J = of the trimer as a function of the atom mass m. To 
guide the eye, dotted lines connect the symbols. 



range (see above) is about 1.5m e , thus providing us with 
an estimate of the error made when extracting the criti- 
cal mass from a set of energies, which excludes the very 
near-threshold regime. 

None of the dimers considered in this subsection sup- 
ports an excited I = state, where I denotes the or- 
bital angular momentum quantum number. Some of the 
dimers do, however, support rotationally excited states 
with I > 0. Even I states are allowed for bosons and 
fermions with opposite spin, and odd I states for spin- 
aligned fermions. Figure|3{c) shows the two-body bound 
state energies for I = (solid line), I = 1 (dashed line) 
and 1 = 2 (dotted line) as a function of the atom mass m. 
The near-threshold behavior of the I > states is, due 
to the presence of the angular momentum barrier of the 
effective potential, distinctly different from that of the 
I = states (for which the angular momentum barrier 
vanishes). The next subsection discusses the energetics 
of the trimer. 



B. N = 3 

Pluses in Fig. 21 show the absolute value of the 
trimer ground state energies, obtained by solving the 
Schrodinger equation using hyperspherical coordinates 
(see Sec. Ill B|l . as a function of the atom mass m. The 
ground state energies £3 extend over nearly three orders 
of magnitude. The most weakly-bound trimer consid- 
ered with m = 5430m e has a ground state energy of 
—4.8 x 10 _4 cm _1 , and the most strongly-bound trimer 
considered with m = 7500m e has a ground state energy 
of -0.27cm- 1 . 

As in the dimer case, we fit the scaled energies yJTn\E^\ 
to the expansion given in Eq. I|13fl with (i < 2) and 
ml 2 ' replaced by cp' and ml 3 ' , respectively. The result- 



ing fitting parameters cp' , Cj 3 ' and ml 3 ' are given in Ta- 
ble||J The fit, shown by a solid line in Fig.[5fb), describes 
the trimer energies quite well. Since the trimer contains 

(3) 

three "dimer bonds" , the critical mass m* for the trimer 
is significantly smaller than that for the dimer. The pa- 
rameter Cq , which quantifies the non-linear dependence 
of the scaled energies on the mass m, is about an order 
of magnitude more negative for the trimer than for the 

(3) 

dimer. The negative value of c 2 suggests that the crit- 

(3) 

ical mass ml predicted by our fit is somewhat smaller 
than the true threshold value, which could be determined 
more precisely if we were able to calculate accurate ener- 
gies closer to threshold. The smallest mass for which we 

reliably determine a negative trimer energy provides an 

(3) 

upper bound for the critical mass m* . We believe that 
the lower bound, i.e., the critical mass predicted by our 
fit, is more accurate than this upper bound. 

Diamonds in Fig. 0] show the first excited state energy 
EP with J = as a function of m. Although the mass at 
which the excited state becomes unbound is larger than 
that at which the ground state becomes unbound, the 
excited state energies approach the threshold in quali- 
tatively the same way as the ground state energies do. 
In fact, the excited state of the helium trimer, which is 
an Efimov state 0, was for a long time considered to 
be the possibly most promising candidate for observing 
universal behaviors experimentally 01 • Just as the near- 
threshold behavior of the excited dimer states with I > 
is qualitatively different from that of the I = state [see 
Fig. Etc)], the near-threshold behavior of trimer states 
with J is predicted to be qualitatively different from 
that of the J = states 0. 



C. N < 40 

We now turn to the discussion of weakly-bound clusters 
with up to N = 40 atoms. Symbols in Fig.^a) show the 
absolute value of the ground state energies Ejy/N per 
particle as a function of the atom mass m for N = 2 — 10. 
The energies for N = 2 and 3 are also shown in Figs. El 
and respectively. The statistical uncertainties of the 
DMC energies Ejy/N per particle, N > 4, are not shown 
in Fig.[3]since they are smaller than the symbol sizes. The 
overall behavior of the ground state energies is similar for 
all N. Below, we use our energies for N = 2 — 10 (see 
Fig. |SJ , and N = 20 and 40 (not shown) to determine 

the critical masses ml^' . 

Symbols in Fig. E^b) show the scaled energies 
\J tti\En\/N for N = 2 — 10 as a function of the atom 



mass m. To determine the critical masses m 



(N) 



we fit 



the scaled energies y/m\EN\ for each TV to the functional 
form given in Eq. (|13fl with c' 2 ' (i < 2) and ml 2 ' re- 
placed by c\ N ^ and to* , respectively. The resulting 
fitting parameters c[ N \ c^' and ml^' are given in Ta- 



7 




FIG. 5: (a) Symbols show the absolute value of the numeri- 
cally determined ground state energies En/N per particle as 
a function of the atom mass m for N = 2 — 10; the energies 
for N = 2 and 3 are also shown in Figs. [3] and 0] respectively. 
To guide the eye, dotted lines connect energies for the same 
N. (b) Symbols show the scaled energies \J m\E^\/N as a 
function of the atom mass for N = 2 — 10. Solid lines show 
our fits to the scaled energies treating c[ , an d m« as 
fitting parameters. For each N, the crossing point of the fit 
with the zero-energy line predicts the critical mass m* . 

ble [U c[ N ^ increases with increasing N, and mi de- 
creases with increasing N. The fitting parameter is 
negative for all AT considered; its largest uncertainty is 
about 40% for N = 20. As discussed in Sec. IfflTfl the 
critical masses m*^ predicted by our fits are expected 
to be lower bounds to the exact threshold value. An up- 
per bound is given by the smallest mass for which we 
report a bound state. A more precise extrapolation of 
the threshold value is complicated by the fact that the 
DMC energies have error bars and that the determina- 
tion of the energy becomes numerically more demanding 
the closer the system's mass is to the critical mass m* . 

Asterisks in Fig. EJa) show the critical masses mi 
predicted by our fits for N = 2 — 10, 20 and 40 as a 
function of 1/N. A diamond in Fig, a) shows the upper 

bound for the critical mass ml°°^ of the bulk system, 
i.e., ml 00 - 1 = 2311. 0TO e (see Sec. Ill All . We choose the 
1/AT-scale since it allows the critical mass for the dimer 
and the bulk system to be shown on the same graph; to 



FIG. 6: (a) Asterisks show the critical masses m» , pre- 
dicted by our fits to the scaled ground state energies \J m|_Ejv|, 
as a function of l/N for TV = 2 — 10, 20 and 40. The dia- 
mond shows an upper bound for the critical mass mi°°' of 
the bulk system and the solid line shows our three-parameter 
fit, Eq. d), with D = 2676(47)m e , E = 11325(458)m e and 
F — — 9696(852)m e to the critical masses mi (the numbers 
in round brackets denote the uncertainty of the fit, neglecting 
possible uncertainties of the critical masses). The dotted line 
shows a lower bound for mi using the equal sign in the re- 
lationship mi N) > mi N ~ 1} (N - 1)/N and mi 2) = 5933.4m e . 
As expected, our numerically determined critical masses lie 
above this analytical bound, (b) Asterisks show the critical 
scattering length ai as a function of 1/N for iV = 2 — 10, 20 
and 40. To guide the eye, a dotted line connects the symbols. 

the best of our knowledge, the functional dependence of 
mi on the system size is unknown. Our critical mass for 
A^ = 40 is significantly larger than the upper bound mf°^ 
determined variationally for the bulk system. Indeed, a 
three-parameter fit of the form 

mi N) =D + E/N + F/N 2 (14) 

to our critical masses for up to 40 atoms, shown by a 
solid line in Fig. EJa), predicts a larger critical mass for 
the bulk system than the upper bound mi°°^ — 2311.0m e 
at which the scattering length crosses zero. We speculate 
that our calculations for comparatively small A^ cannot 
be used to extrapolate mi 00 ^ reliably since the ratio of 
bulk to surface atoms increases appreciably with increas- 
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ing N. Furthermore, our extrapolated critical masses 
have non-negligible uncertainties. To connect the results 
discussed here to a realistic physical system, we note that 
homogeneous atomic spin-polarized hydrogen, which has 
an atom mass of m = 1837m e and interacts through a 
sum of two-body potentials only slightly different from 
that considered here, exists under normal pressure as a 
gas and not as a liquid 45] . 

The critical mass mi N ^ of the TV-body system is 
bounded by the critical mass of the system with N — 1 
particles through mi N) > mi N ^ (N - 1)/N [I3|. A dot- 
ted line in Fig. Ela) shows this lower bound, assuming 

mi^ = 5933. 4m e . Figure Etb) indicates that this ana- 
lytical estimate provides a weak bound for all bosonic 
systems considered here. An upper bound is given by 

(TV) (JV-1) 

mi — mi 

To relate our critical masses to an experimentally tun- 
able parameter, we calculate the scattering length for 
each mi N ^ and refer to it as the critical scattering length 
a* . Asterisks in Fig. Efb) show the critical scattering 
length a[ NS> as a function of 1/N. Figure HJb) suggests 
that Borromcan trimers exist for a < — 58.8ao and Bor- 
romean tetramers for a £ [— 58.8ao, — 26.0ao]. Investiga- 
tion of the stability of these weakly-bound Borromean 
states is beyond the scope of this paper. 

D. Correlations 

We now investigate correlations between energies of 
the three- and four-particle systems. The description of 
universal properties of the trimer Q, such as the de- 
scription of Efimov states, requires two parameters, a 
two-body momentum scale (typically taken to be in- 
versely proportional to the s-wave scattering length a) 
and a three-body momentum scale /i^ 3 ' (in some studies, 
the three-body parameter A* is used instead 0> 
While the universal behaviors of the trimer are quite 
well understood 9], much less is known about those of 
larger systems. For example, although one expects a new 
momentum scale fj,^ to be needed for the description 
of universal properties of the tetramer 0, 0] , there is 
evidence that at least some observables of the tetramer 
near threshold are independent of this new momentum 
scale 0. In this context, a number of studies have fo- 
cused on the Tjon line 0, El H3, El E2 , which was first 
investigated in nuclear physics and refers to the approxi- 
mately linear correlation between the energies of the four- 
nucleon and the three-nucleon system |5ll |52j. In the 
following, correlations between our trimer and tetramer 
energies, which are calculated for realistic atom-atom in- 
teractions, are demonstrated. 

Pluses in Figs. Efa) and (b) show the ratio between 
the ground state energies E4 and E2 as a function of the 
ratio between the ground state energies E3 and E2 on a 
linear and double-logarithmic scale, respectively. Unlike 
in Tjon's original work for fixed dimer energy [5lL l52| . 




10° 10 1 10 2 10 3 10 4 

E 3 /E 2 

FIG. 7: Pluses show the energy ratio E4/E2 as a function of 
the energy ratio E3/E2 (E4, E3 and E2 denote ground state 
energies) on (a) a linear scale and (b) a log-log scale. Solid 
lines show our two-parameter fit (see text). Dotted lines show 
the result obtained within an effective quantum mechanical 
approach The data with the largest energy ratios cor- 
respond to systems closest to threshold, and those with the 
smallest energy ratios correspond to systems farthest away 
from threshold. 



we scale the trimer and tetramer energies in Fig. [7| by 
the dimer energy since E2 depends on the atom mass m. 
For each data point, the ground state energies E2, E3 
and E4 are calculated for the same atom mass m. The 
systems closest to threshold are those with the largest 
energy ratios. For the smallest mass, m = 5950m e , in- 
cluded in Fig.[7Jb), the absolute value of E2 is nearly four 
orders of magnitude smaller than that of E3, and more 
than four orders of magnitude smaller than that of E4. 
A two-parameter fit of the form E4/E2 = B 3 + C 3 E 3 /E2, 
shown by solid lines in Figs.[7{a) and (b), describes the 
dependence of E4/E2 on E 3 /E 2 quite well (especially for 
systems close to threshold) , thus confirming the existence 
of the Tjon line for atomic clusters. In particular, we find 
B 3 = -34.9(8.3) and C 3 = 5.008(5) (see also Table HQ, 
where the numbers in brackets indicate the uncertainty 
of the fit, neglecting possible inaccuracies of the energy 
ratios. 

For comparison, dotted lines in Figs. Eta) and (b) show 
a result derived within an effective quantum mechanics 
approach applied to bosonic clusters with N = 2 — 4 
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helium atoms gj. This study finds B 3 — -24.752 and 
C 3 = 4.075 for 69 < E 3 /E 2 < 142. This range of E 3 /E 2 
values is significantly smaller than that considered in the 
present paper [solid lines in Figs. E{a) and (b)]. The 
slope of the Tjon line derived within the effective quan- 
tum mechanical approach, applied to helium clusters, is 
somewhat smaller than our slope, which is derived from 
a series of numerical calculations for weakly-bound Van 
der Waals clusters. We find that our slope decreases if 
we perform a fit that excludes energy ratios for systems 
very close to threshold. This may explain the discrepancy 
between the results obtained within the two approaches. 

It has been argued that, if the four-body momen- 
tum scale fi^' coincides with the three-body momen- 
tum scale fi^ 14], the slope of the Tjon line is about 
five. This argument suggests that the systems studied in 
the present paper have approximately equal three- and 
four-body momentum scales. It has further been sug- 
gested that the existence of equal three- and four-body 
momentum scales can be traced back to the repulsive 
core of the two-body interactions 0] . This interpreta- 
tion suggests that the Tjon line is roughly five for all 
atomic systems near threshold and that the ground state 
energy of any weakly-bound tetramer interacting addi- 
tively through Van der Waals potentials with repulsive 
core can be estimated quite reliably if the corresponding 
dimer and trimer ground state energies are known. This 
interpretation could be checked more rigorously by per- 
forming a series of calculations for systems interacting ad- 
ditively through a shape-dependent two-body potential, 
which depends on a parameter that controls the softness 
of the repulsive short-range part of the potential. For 
a "soft" repulsive core, the four-body momentum scale 
should deviate from the three-body scale and the slope 
of the Tjon line should deviate from five. Such a study 
is beyond the scope of this paper. 

We now consider correlations between the tetramer 
ground state energy E± and the trimer excited state en- 
ergy E$ , both scaled by the dimer ground state en- 
ergy E 2 . Pluses in Fig. [S] show the energy ratio E4/E2 
as a function of the energy ratio E% /E<z on a linear 
scale. A solid line shows our two-parameter fit E4/E2 = 
B 3 + C 3 E [ p /E 2 with B 3 = -558(35) and C 3 = 565(29), 
while a dotted line shows that derived in Ref. for 
1.54 < E^ /E 2 < 2 with B 3 = -742.0 and C 3 = 645.1. 
The agreement of the slopes, derived within two differ- 
ent frameworks and applied to two different systems, is 
quite reasonable. We now use the linear dependence of 
E 4 /E 2 on E 3 /E 2 and of E A /E 2 on E^ /E 2 to predict the 
slope C3 for the linear dependence of E 3 /E 2 on E% / 'E 2 , 
E 3 /E 2 = B 3 + C 3 E i 3 1) /E 2 . From our slopes C 3 and C a , 
we obtain C 3 — 112(6). For comparison, a fit to our data 
gives C 3 = 124(5). The good agreement lends support 
to the predictive power of the approximately linear cor- 
relations between energy ratios of clusters with varying 
number of atoms. 
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FIG. 8: Pluses show the energy ratio E4/E2 as a function of 
the energy ratio / E2 on a linear scale (E2 and E4 denote 
ground state energies, and denotes the first excited state 
energy with J = 0). A solid line shows our two-parameter fit 
(see text). A dotted line shows the results obtained within an 
effective quantum mechanical approach ||. The data with the 
largest energy ratios correspond to systems closest to thresh- 
old, and those with the smallest energy ratios correspond to 
systems farthest away from threshold. 

Correlations between the energies of two clusters dif- 
fering in size by one atom, i.e., a linear relationship of 
the form Ej^+i = bjy + cnEn, have been predicted an- 
alytically based on a separable approximation scheme 
for any cluster size |53j | and numerically by performing 
variational calculations for small mixed 3 Hei- 4 He., clus- 
ters with i + j < 5 [H(J . We find that our energies of 
clusters differing in size by one atom are not well de- 
scribed by such a linear two-parameter fit. If we instead 
scale, as in the investigation of the correlations between 
the tetramer and trimer energies, our ground state en- 
ergies of the (N + 1)- and V-atom clusters by the en- 
ergy of the (N — l)-atom cluster, we find an approxi- 
mately linear relationship. Pluses in Figs. EJa)-(f) show 
the energy ratios E^+i/ '-Ejv-i as a function of En / ' E^-i 
for N = 4 — 9, and solid lines a linear fit of the form 
En+i/En-i = Bn + CnEn/En-i. The fitting parame- 
ters Bn and Cm are summarized in Table ITT1 We refer to 
the approximately linear dependence of the energy ratios 
of clusters, which is illustrated in Fig. ^\ as generalized 
Tjon lines. It will be interesting to investigate the impli- 
cations of the behaviors of these generalized Tjon lines 
for the universal properties of weakly-bound bosonic clus- 
ters. 



IV. STRUCTURAL PROPERTIES 

This section presents selected structural properties of 
weakly-bound bosonic clusters in their ground state. The 
expectation values for N > 4 are calculated using the MC 
estimator given in Eq. (|10fl . Figure shows the pair 
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FIG. 9: Pluses show the energy ratio Em+\/ En-\ as a func- 
tion of the energy ratio En/En-i for (a) N = 4, (b) N = 5, 
(c) N = 6, (d) N = 7, (e) TV = 8, and (f) N = 9. Solid lines 
show linear fits of the form En+i/En-i = Bn+CnEn/En-i 
(see text). Note that the range of the vertical axis extends 
from 1 to 8 in all panels while that of the horizontal axis 
varies. 



distribution function P(r) for N = 4 and four different 
masses, i.e., m = 5950m e , 5750m e , 5400m e and 5150m e . 
The pair distribution function P(r) indicates the likeli- 
hood of finding two particles at a distance r from each 
other and is normalized so that 



P(r)r 2 dr = 1. 



(15) 



As the mass decreases, the maximum of P{r), which is 
located at r « 10ao, decreases. Furthermore, the pair 
distribution functions extend to significantly larger r val- 
ues for small m than for large m. For example, the 
largest interparticle distance sampled in our DMC runs 
for m = 5950m e is r w lOOao while that for m — 5150m e 
is r w 200ao- We find that the densities, not shown, 
of the weakly-bound clusters studied in this paper are 
structureless and do not possess any shell structure. The 
highly-diffuse clusters can thus be most appropriately 
thought of as "diffuse liquid blobs" . 

To compare the structural behaviors of clusters with 
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73 jv 


Cjv 


3 


-34.9(8.3) 5.008(5) 


4 


-3.37(27) 


3.10(8) 


5 


-2.27(08) 


2.78(2) 


6 


-1.71(03) 


2.49(2) 


7 


-1.55(07) 


2.45(2) 


8 


-1.50(06) 


2.36(3) 


9 


-1.11(11) 


2.18(5) 



TABLE II: Fitting parameters B N and Cjv for N = 3 - 
9. Numbers in brackets denote the uncertainty of the two- 
parameter fit En+\/ En-i = Bn + CnEn/En-i, neglecting 
possible uncertainties of the energy ratios. 
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FIG. 10: Pair distribution functions P(r) as a function of the 
interparticle distance r for N = 4 and four different masses, 
i.e., m — 5950m e , 5750m e , 5400m e and 5150m e (from top to 
bottom). The pair distribution functions are calculated using 
the MC estimator given in Eq. (| 1 0|l , which combines the VMC 
and DMC expectation values. Statistical uncertainties are not 
shown for clarity. 



different TV, symbols in Fig. ^] show the expectation 
values of the interparticle distance r for clusters in the 
ground state with N = 2—10 (denoted by (r)jv) as a func- 
tion of the atom mass m. For fixed N, (r) n decreases, 
as expected, with decreasing mass m. For a given mass, 
(r)jv decreases with increasing N. This behavior is con- 
sistent with the fact that the energy per particle for fixed 
mass decreases with increasing N. Furthermore, for fixed 
m, the expectation values (r) jy should reach a constant 
in the large N limit. Indeed, Fig. ^] indicates that the 
difference between (r) n for two clusters differing in size 
by one atom is smaller for large than for small N. 

It has been suggested that scaling functions, which al- 
low the structural properties of the tetramer to be ex- 
pressed in terms of expectation values of the dimer and 
trimer, exist |14| : the exact functional forms are, how- 
ever, to the best of our knowledge unknown. Symbols 
in Fig. 1121 show the ratio (r)i/(r)2 as a function of the 
ratio (r}s/ (r)2- For each data point, the expectation val- 
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gases. Indeed, the observation of resonances in an ultra- 
cold Bose gas has recently been interpreted as evidence 
for the presence of loosely-bound Efimov trimers [2l| . 
Resonances associated with tetramers have also been re- 
ported j22j. These experiments may just be the begin- 
ning of detailed studies of the rich behaviors of few-body 
systems under controlled conditions. On the theoretical 
side, little is known about the universal near-threshold 
behaviors of three-dimensional systems with more than 
three particles. The reason is that analytical treatments 
become increasingly more complex as the number of de- 
grees of freedom increases. On the other hand, numerical 
treatments are complicated by the fact that the kinetic 
and potential energy nearly cancel, thus requiring both 
to be calculated with high accuracy. 



FIG. 11: Expectation value (r) n of the interparticle distance 
for clusters in the ground state with TV = 2 — 10 (from right to 
left) as a function of the atom mass m. The expectation values 
for TV > 4 are calculated using the extrapolated estimator, 
Eq. ill (H . which combines the VMC and DMC expectation 
values. Errorbars, not shown, are at most about three times 
as large as the symbol sizes. 
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<r> 3 /<r> 2 

FIG. 12: Ratio {r)4,/{r)2 as a function of the ratio (r)z/{r)2- 
Pluses show our numerically determined ratios while the solid 
line shows a two-parameter fit (see text). 

ues (r)i, (r)a and {r) 2 are calculated for the same mass. 
The errorbars, not shown, are smaller than twice the size 
of the symbols. These ratios are well described by a 
two-parameter fit of the form (r)i/(r)2 = G(r)a/(r)2 + 
H((r) 3 /(r) 2 ) 2 with G = 0.751(5) and H = 0.313(7) (solid 
line). We hope that the structural properties presented 
here will stimulate and aid further studies of weakly- 
bound bosonic clusters. 



This paper presents a detailed study of the near- 
threshold behaviors of weakly-bound three-dimensional 
bosonic clusters with up to 40 atoms, for which the un- 
derlying potential energy surface is written as a sum of 
realistic Van der Waals atom-atom potentials with short- 
range repulsion and attractive long-range tail. In particu- 
lar, we determine the critical mass m* for clusters with 
/V = 2 — 10, 20 and 40 by performing calculations for each 
cluster as a function of the atom mass m. To the best 
of our knowledge, these are the first calculations that at- 
tempt an accurate determination of the critical coupling 
strengths of clusters with up to 40 atoms. Our critical 
masses are compared to analytical bounds. Furthermore, 
we show that our numerically determined three- and four- 
particle energies, scaled by the corresponding dimer en- 
ergies, fall on the Tjon line. We present numerical evi- 
dence that the scaled energies of larger clusters differing 
in size by one atom also correlate linearly, i.e., the en- 
ergy ratios fall on what we refer to as generalized Tjon 
lines. Motivated by a recent calculation based on zero- 
range models 0], we speculate that all atomic cluster 
systems show similar near-threshold behaviors. Finally, 
we present selected structural properties of weakly-bound 
few-body systems. 



In closing, we emphasize that the near-threshold be- 
havior of clusters crucially depends on the dimensional- 
ity. For example, the near-threshold behavior of weakly- 
bound two-dimensional few-body systems 54, 55] is very 
different from that presented here for three-dimensional 
systems. We hope that our work will stimulate further 
experimental and theoretical work on weakly-bound clus- 
ters. 



V. CONCLUSION 

The physics of weakly-bound few-body systems can 
experimentally be investigated using ultracold atomic 
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